home *** CD-ROM | disk | FTP | other *** search
/ Space & Astronomy / Space and Astronomy (October 1993).iso / mac / programs / space / DE118I.ZIP / SSYSTEM.C < prev    next >
C/C++ Source or Header  |  1993-03-01  |  17KB  |  834 lines
  1. /* Main program for numerical integration of the solar system.
  2.  * See the readme file for discussion.
  3.  *
  4.  * Before compiling, set type of computer in mconf.h,
  5.  * arithmetic precision in prec.h, and model parameters
  6.  * in ssystem.h.
  7.  *
  8.  * If DE200 compatible outputs are desired, include the
  9.  * rotation r118_200() in the output routine.
  10.  */
  11.  
  12. /* install trap handler for arithmetic debugging */
  13. #define FPESHOW 0
  14.  
  15. #ifndef NOINCS
  16. #include "mconf.h"
  17. #include "prec.h"
  18. #include "ssystem.h"
  19. #endif
  20. #include "int.h"
  21. #include "const.h"
  22.  
  23. #ifdef DEBUG
  24. #undef DEBUG
  25. #endif
  26. #define DEBUG 0
  27.  
  28. extern DOUBLE GMs[], KG, C, EMRAT;
  29. extern DOUBLE yn0[], JD0, yn1[], JD1;
  30.  
  31. /* Define 1 for ephemeris output in double precision even if LDOUBLE=1.
  32.  * The saved integrator state will still be LDOUBLE.
  33.  */
  34. #define OUTDOUBLE 1
  35.  
  36. DOUBLE vnewt[6*NTOTAL];
  37. DOUBLE *awork; /* Adams integrator work array */
  38. DOUBLE *rwork; /* Runge-Kutta work array */
  39. DOUBLE Rij[NTOTAL][NTOTAL];
  40. DOUBLE Rij3[NTOTAL][NTOTAL];
  41.  
  42. int fd;
  43. char fname[80];
  44.  
  45. #if MSC
  46. #include <stdio.h>
  47. #include <sys\types.h>
  48. #include <sys\stat.h>
  49. #include <fcntl.h>
  50. #endif
  51. #if UNIX
  52. #include <stdio.h>
  53. #include <sys/types.h>
  54. #include <sys/stat.h>
  55. #include <fcntl.h>
  56. #endif
  57. #if VAX
  58. #include <stdio.h>
  59. #include <file.h>
  60. #include <stat.h>
  61. #endif
  62. #if RSX
  63. #include <stdio.h>
  64. #include <fcntl.h>
  65. #define CPERM 1
  66. #endif
  67. #if SVC
  68. #endif
  69.  
  70.  
  71. DOUBLE ssyt, err, ssyh, unused;
  72. int ord, outsteps;
  73. long isamps, nsamps, nsteps;
  74.  
  75. /* Subroutines to save and restore the integrator state
  76.  */
  77. void savstate()
  78. {
  79. /*int f, i;*/
  80. register int f;
  81.  
  82. #if MSC
  83. f = open( "de118.sav", O_BINARY | O_CREAT | O_RDWR, S_IREAD | S_IWRITE );
  84. #endif
  85. #if UNIX
  86. f = open( "de118.sav", O_CREAT | O_TRUNC | O_RDWR, S_IWRITE | S_IREAD );
  87. #endif
  88. #if VAX
  89. f = open( "de118v.sav", O_CREAT | O_TRUNC | O_WRONLY, S_IREAD | S_IWRITE, "rfm = var" );
  90. #endif
  91. #if RSX
  92. f = open( "de118.sav", O_CREAT | O_BIN | O_RDWR, CPERM );
  93. #endif
  94. #if SVC
  95. f = creat( "de118.sav" );
  96. #endif
  97. sinternal(f);
  98. write( f, &ssyt, sizeof(DOUBLE) );
  99. write( f, &err, sizeof(DOUBLE) );
  100. write( f, &ssyh, sizeof(DOUBLE) );
  101. write( f, &unused, sizeof(DOUBLE) );
  102. write( f, yn0, 6*NTOTAL*sizeof(DOUBLE) );
  103. write( f, &ord, sizeof(int) );
  104. write( f, &outsteps, sizeof(int) );
  105. write( f, awork, ((NEQ+1)*(MAXORD+2)+(6*NEQ))*sizeof(DOUBLE) );
  106. /*write( f, awork, i );*/
  107. close(f);
  108. }
  109.  
  110. void resstate()
  111. {
  112. /*int f, i;*/
  113. register int f;
  114.  
  115. #if MSC
  116.     f = open( "de118.sav", O_BINARY | O_RDONLY, S_IREAD );
  117. #endif
  118. #if UNIX
  119.     f = open( "de118.sav", O_RDONLY, S_IREAD | S_IWRITE );
  120. #endif
  121. #if VAX
  122.     f = open( "de118.sav", O_RDONLY, S_IREAD );
  123. #endif
  124. #if RSX
  125.     f = open( "de118.sav", O_BIN | O_RDONLY, CPERM );
  126. #endif
  127. #if SVC
  128.     f = open( "de118.sav", 0 );
  129. #endif
  130. if( f <= 0 )
  131.     {
  132.     printf( "Can't find de118.sav\n" );
  133.     exit(0);
  134.     }
  135. rinternal(f, awork);
  136. read( f, &ssyt, sizeof(DOUBLE) );
  137. read( f, &err, sizeof(DOUBLE) );
  138. read( f, &ssyh, sizeof(DOUBLE) );
  139. read( f, &unused, sizeof(DOUBLE) );
  140. read( f, yn0, 6*NTOTAL*sizeof(DOUBLE) );
  141. read( f, &ord, sizeof(int) );
  142. read( f, &outsteps, sizeof(int) );
  143. read( f, awork, ((NEQ+1)*(MAXORD+2)+(6*NEQ))*sizeof(DOUBLE) );
  144. /*read( f, awork, i );*/
  145. close(f);
  146. }
  147.  
  148. /* Double precision variables for printf() */
  149. static double dx, dy, dz;
  150.  
  151. main()
  152. {
  153. DOUBLE xa, ya, za, xb, yb, zb; /* for display */
  154. DOUBLE *pdv;
  155. int i, ii, iout, restored;
  156. DOUBLE adstep();
  157.  
  158. #if FPESHOW
  159. fpeins();
  160. #endif
  161. ii = (NEQ+1)*(MAXORD+2)+(7*NEQ);
  162. awork = (DOUBLE *)malloc( ii * sizeof(DOUBLE) );
  163. if( awork == 0 )
  164.     {
  165.     printf( "malloc(awork) failed\n" );
  166.     exit(0);
  167.     }
  168. /* Note, it is not really necessary to clear out the work arrays.
  169.  */
  170. pdv = awork;
  171. xa = Zero;
  172. for( i=0; i<ii; i++ )
  173.     *pdv++ = xa;
  174. printf( "Allocated %d bytes to awork\n", i );
  175.  
  176. i = (MAXRUNG+3)*NEQ*sizeof(DOUBLE);
  177. rwork = (DOUBLE *)malloc(i);
  178. if( rwork == 0 )
  179.     {
  180.     printf( "malloc(rwork) failed\n" );
  181.     exit(0);
  182.     }
  183. printf( "Allocated %d bytes to rwork\n", i );
  184.  
  185. printf( "NEQ %d, NEQC %d, 6*NTOTAL %d, IAROIDS %d, NBODY %d, NTOTAL %d\n",
  186.    NEQ, NEQC, 6*NTOTAL, IAROIDS, NBODY, NTOTAL );
  187.  
  188. /* Read in all physical constants from file named aconst.h.
  189.  * This is optional, since the constants are already compiled
  190.  * in to the program.
  191.  */
  192. rdnums();
  193.  
  194. restored = 0;
  195. printf( "Do you want to restore a previously saved integrator state?" );
  196. /* If the answer is 'z' then the state will be restored and also
  197.  * the program will ask for a new interval between output records.
  198.  */
  199. scanf( "%c", fname );
  200. printf( "%s\n", fname );
  201. ii = fname[0] & 0x7f;
  202. if( (ii == 'y') || (ii == 'z') )
  203.     {
  204.     resstate();
  205.     restored = 1;
  206.     goto opnout;
  207.     }
  208.  
  209. orlup:
  210. printf( "Adams order ? " );
  211. scanf( "%d", &ord );
  212. printf( "%d\n", ord );
  213. if( (ord > MAXORD) || (ord < 1) )
  214.     {
  215.     printf( "order must be between 1 and %d\n", MAXORD );
  216.     goto orlup;
  217.     }
  218.  
  219. printf( "step size, days ? " );
  220. scanf( "%lf", &dx );
  221. printf( "%.8lf\n", dx );
  222. ssyh = dx;
  223.  
  224.  
  225. printf( "initializing...\n" );
  226.  
  227. /* Ensure that the relativistic center of the N bodies is at 0.0.
  228.  * This condition is already approximately satisfied by states
  229.  * taken from the JPL ephemeris tapes.
  230.  */
  231. findcenter( JD0, yn0 );
  232. findcenter( JD1, yn1 );
  233.  
  234. ssyt = JD0;
  235. err = Zero;
  236. rkstart( NEQ, rwork );
  237. #if DEBUG
  238. printf( "rkstart " );
  239. #endif
  240. adstart( ssyh, &yn0[0], awork, NEQ, ord, ssyt );
  241. #if DEBUG
  242. printf( "adstart\n" );
  243. #endif
  244. printf( "initialized.\n" );
  245.  
  246.  
  247.  
  248. /* Each output file record contains the Julian date
  249.  * followed by the current velocity and position states.
  250.  * See ini118.h for the data structure of the state array.
  251.  * The state variables for a test object such as a comet
  252.  * or asteroid should be placed in the array immediately after
  253.  * the Sun state.
  254.  */
  255. opnout:
  256.  
  257. /* If a previous state was restored, you had the option of 
  258.  * getting the interval between outputs from the saved state
  259.  * file or of entering a new interval by responding 'z'.
  260.  */
  261. if( ii != 'y' )
  262.     {
  263.     printf( "Interval between output samples, in steps ? " );
  264.     scanf( "%d", &outsteps );
  265.     printf( "%d\n", outsteps );
  266.     }
  267.  
  268. printf( "Output file name ? " );
  269. scanf( "%s", fname );
  270. printf( "%s\n", fname );
  271. #if MSC
  272. fd = open( fname, O_BINARY | O_CREAT | O_RDWR, S_IREAD | S_IWRITE );
  273. #endif
  274. #if UNIX
  275. fd = open( fname, O_CREAT | O_TRUNC | O_RDWR, S_IWRITE | S_IREAD );
  276. #endif
  277. #if VAX
  278. fd = open( fname, O_CREAT | O_TRUNC | O_WRONLY, S_IREAD | S_IWRITE, "rfm = var" );
  279. #endif
  280. #if RSX
  281. fd = open( fname, O_CREAT | O_BIN | O_RDWR, CPERM );
  282. #endif
  283. #if SVC
  284. fd = creat( fname );
  285. #endif
  286. if( fd <= 0 )
  287.     {
  288.     printf( "Can't create \"%s\"\n", fname );
  289.     exit(1);
  290.     }
  291.  
  292. printf( "Number of output samples ? " );
  293. scanf( "%ld", &nsamps );
  294. printf( "%ld\n", nsamps );
  295.  
  296. /* Write initial state vector to output file */
  297. if( restored == 0 )
  298.     wfile( ssyt, yn0 );
  299.  
  300.  
  301. nsteps = 0;
  302. for( isamps=0; isamps<nsamps; isamps++ )
  303.     {
  304.     for( iout=0; iout<outsteps; iout++ )
  305.         {
  306.         err += adstep( &ssyt, &yn0[0], NEQC );
  307.         nsteps += 1;
  308.         }
  309.     wfile( ssyt, yn0 );
  310.     dx = err;
  311.     dy = ssyt;
  312.     printf( "%5ld %11.2lf %.6e\n", nsteps, dy, dx );
  313. #if IBMPC
  314. /* Check for operator abort
  315.  */
  316.     if( kbhit() )
  317.         {
  318.         ii = getchar();
  319.         if( (ii & 0x7f) == 'S' )
  320.             {
  321.             printf( "Operator interrupt, closing %s\n", fname );
  322.             break;
  323.             }
  324.         }
  325. #endif
  326.     }
  327. close( fd );
  328. savstate();
  329.  
  330. /* Test if the ending time is equal to the time of
  331.  * the compiled-in test state vector. If so, print out
  332.  * the errors.
  333.  */
  334. xa = ssyt - JD1;
  335. if( xa < 0 )
  336.     xa = -xa;
  337.  /* allow for roundoff in the time variable */
  338. if( xa > 1e-6 )
  339.     iout = 0; /* error printout not appropriate */
  340. else
  341.     iout = 1;
  342.  
  343. printf( "Final x, y, z, positions and errors:\n" );
  344. /*
  345. ii = 6*FMASS;
  346. for(i=FMASS; i<IAROIDS; i++)
  347. */
  348.  
  349. ii = 0;
  350. for(i=0; i<=ISUN; i++)
  351.  
  352.     {
  353. /* Compare Heliocentric positions
  354.     xa = yn0[ii+1] - yn0[6*ISUN+1];
  355.     ya = yn0[ii+3] - yn0[6*ISUN+3];
  356.     za = yn0[ii+5] - yn0[6*ISUN+5];
  357.     xb = yn1[ii+1] - yn1[6*ISUN+1];
  358.     yb = yn1[ii+3] - yn1[6*ISUN+3];
  359.     zb = yn1[ii+5] - yn1[6*ISUN+5];
  360. */
  361. if( i == 0 )
  362.     {
  363.     dx = yn0[0];
  364.     dy = yn0[2];
  365.     dz = yn0[4];
  366.     printf("%19.15lf%19.15lf%19.15lf\n", dx, dy, dz );
  367.     }
  368.  
  369.     xa = yn0[ii+1] - yn1[ii+1];
  370.     ya = yn0[ii+3] - yn1[ii+3];
  371.     za = yn0[ii+5] - yn1[ii+5];
  372.     err = SQRT( xa*xa + ya*ya + za*za );
  373. /* Use ieee.c to print numerical results */
  374. #if LDOUBLE
  375.     e64toasc( &yn0[ii+1], fname, 15 );
  376.     printf( "%s ", fname );
  377.     e64toasc( &yn0[ii+3], fname, 15 );
  378.     printf( "%s ", fname );
  379.     e64toasc( &yn0[ii+5], fname, 15 );
  380.     printf( "%s ", fname );
  381.     if( iout )
  382.         {
  383.         e64toasc( &err, fname, 3 );
  384.         printf( "%s\n", fname );
  385.         }
  386. #else
  387.     e53toasc( &yn0[ii+1], fname, 15 );
  388.     printf( "%s ", fname );
  389.     e53toasc( &yn0[ii+3], fname, 15 );
  390.     printf( "%s ", fname );
  391.     e53toasc( &yn0[ii+5], fname, 15 );
  392.     printf( "%s ", fname );
  393.     if( iout )
  394.         {
  395.         e53toasc( &err, fname, 3 );
  396.         printf( "%s\n", fname );
  397.         }
  398. #endif
  399. /*
  400.     dx = (double) yn0[ii+1];
  401.     dy = (double) yn0[ii+3];
  402.     dz = (double) yn0[ii+5];
  403.     printf("%19.15lf%19.15lf%19.15lf", dx, dy, dz );
  404.     dz = (double) err;
  405.     printf(" %.3e\n", dz );
  406. */
  407.     ii += 6;
  408.     }
  409. #if FPESHOW
  410. fperem();
  411. #endif
  412. } /* end of main program */
  413.  
  414.  
  415.  
  416. /* Subroutine func() calculates the forces and accelerations.
  417.  * Data in the output vector v[] are in the order
  418.  * d^2x/dt^2, dx/dt, d^2y/dt^2, dy/dt, d^2z/dt^2, dz/dt
  419.  * for each object.  For this problem the velocities dx/dt, ...
  420.  * are simply copied over from the input array yw[].
  421.  */
  422.  
  423. #if MOON
  424. DOUBLE yw[6*NTOTAL];
  425. DOUBLE ymoon[6];
  426. #endif
  427.  
  428. func( t, yin, v )
  429. DOUBLE t; /* time */
  430. DOUBLE yin[]; /* input state: velocity and position */
  431. DOUBLE v[]; /* output: calculated acceleration, copy of input velocity */
  432. {
  433. DOUBLE xs, ys, zs, xd, yd, zd, xv, yv, zv, temp;
  434. int i, j, ii, jj;
  435.  
  436. #if DEBUG
  437. printf( "func " );
  438. #endif
  439. #if MOON
  440. /* Locally replace input variable EMB by barycentric earth
  441.  * and input variable Moon-Earth by barycentric Moon.
  442.  */
  443. fromemb( yin, yw );
  444. #if DEBUG
  445. printf( "fromemb " );
  446. #endif
  447. #endif
  448.  
  449. /* asteroid positions */
  450. #if AROIDS
  451. aroids( t, yw );
  452. #if DEBUG
  453. printf( "aroids " );
  454. #endif
  455. #endif
  456.  
  457.  
  458. /* Table of distances between objects i and j */
  459. distances(yw);
  460. #if DEBUG
  461. printf( "distances " );
  462. #endif
  463.  
  464. fixsun( yw, v );
  465. #if DEBUG
  466. printf( "fixsun " );
  467. #endif
  468.  
  469. #if MOON
  470. for( j=0; j<6; j++ )
  471.     yin[(6*ISUN)+j] = yw[(6*ISUN)+j];
  472. #endif
  473.  
  474. /* Compute Newtonian acceleration. */
  475. ii = 6*FMASS;
  476. for( i=FMASS; i<NTOTAL; i++ )
  477.     {
  478.     xs = Zero;
  479.     ys = Zero;
  480.     zs = Zero;
  481.     xv = yw[ii+1];
  482.     yv = yw[ii+3];
  483.     zv = yw[ii+5];
  484.     jj = 6*FMASS;
  485.     for( j=FMASS; j<NTOTAL; j++ )
  486.         {
  487.         if( (j != i) && (j != ISUN) )
  488.             {
  489.             xd = yw[jj+1] - xv;
  490.             yd = yw[jj+3] - yv;
  491.             zd = yw[jj+5] - zv;
  492.             temp = GMs[j] * Rij3[i][j];
  493.             xs += xd * temp;
  494.             ys += yd * temp;
  495.             zs += zd * temp;
  496.             }
  497.         jj += 6;
  498.         }
  499. /* Add the Sun last. */
  500.     if( i != ISUN )
  501.         {
  502.         xd = yw[6*ISUN+1] - xv;
  503.         yd = yw[6*ISUN+3] - yv;
  504.         zd = yw[6*ISUN+5] - zv;
  505.         temp = GMs[ISUN] * Rij3[i][ISUN];
  506.         xs += xd * temp;
  507.         ys += yd * temp;
  508.         zs += zd * temp;
  509.         }
  510.     vnewt[ii] = xs; /* acceleration */
  511.     vnewt[ii+2] = ys;
  512.     vnewt[ii+4] = zs;
  513.     vnewt[ii+1] = yw[ii]; /* velocity */
  514.     vnewt[ii+3] = yw[ii+2];
  515.     vnewt[ii+5] = yw[ii+4];
  516.     ii += 6;
  517.     }
  518.  
  519. #if DOREL
  520. reltiv( yw, v );
  521. #if DEBUG
  522. printf( "reltiv " );
  523. #endif
  524. #else
  525. /* No relativity theory. Return the Newtonian accelerations. */
  526. ii = 6*FMASS;
  527. for( i=FMASS; i<IAROIDS; i++ )
  528.     {
  529.     v[ii] = vnewt[ii];
  530.     v[ii+2] = vnewt[ii+2];
  531.     v[ii+4] = vnewt[ii+4];
  532.     v[ii+1] = yw[ii];
  533.     v[ii+3] = yw[ii+2];
  534.     v[ii+5] = yw[ii+4];
  535.     ii += 6;
  536.     }
  537. #endif /* DOREL */
  538.  
  539.  
  540. #if MOON
  541. oblate( t, yw, v, ymoon );
  542. #if DEBUG
  543. printf( "oblate\n" );
  544. #endif
  545.  
  546. /* Add the Newtonian accelerations last. */
  547. ii = 6*FMASS;
  548. for( i=FMASS; i<IAROIDS; i++ )
  549.     {
  550.     v[ii] += vnewt[ii];
  551.     v[ii+2] += vnewt[ii+2];
  552.     v[ii+4] += vnewt[ii+4];
  553. /*
  554.  *    v[ii+1] = yw[ii];
  555.  *    v[ii+3] = yw[ii+2];
  556.  *    v[ii+5] = yw[ii+4];
  557.  */
  558.     ii += 6;
  559.     }
  560. /* Convert barycentric Earth and Moon to output EMB and M variables. */
  561. ii = 6*IEARTH;
  562. jj = 6*IMOON;
  563. for( i=0; i<6; i += 2 )
  564.     {
  565.     xd = v[ii+i]; /* Earth */
  566.     yd = v[jj+i]; /* Moon */
  567.     v[ii+i] = (EMRAT * xd + yd)/(EMRAT+One); /* EMB */
  568.     v[jj+i] = yd - xd; /* M = Moon - Earth */
  569.     v[ii+i+1] = yin[ii+i]; /* copy original velocity */
  570.     v[jj+i+1] = yin[jj+i];
  571.     }
  572. #if 0
  573. /* Display initial acceleration discrepancies (see findcent.c) */
  574. chkacc(v);
  575. exit(0);
  576. #endif
  577.  
  578. #endif /* MOON */
  579. }
  580.  
  581.  
  582. /* Constrain the barycenter to stay at the origin.
  583.  * This is done by adjusting the Sun.
  584.  */
  585.  
  586. DOUBLE mustar[NTOTAL];
  587.  
  588. fixsun( y, v )
  589. DOUBLE y[], v[];
  590. {
  591. DOUBLE xx, yy, zz, vx, vy, vz;
  592. DOUBLE csqi, s;
  593. int i, j, k, ii, jj;
  594.  
  595.  
  596. csqi = Half / (C*C);
  597. for( k=0; k<2; k++ )
  598.  { /* Iterate to find solution of implicit expressions. */
  599.  
  600.  
  601. /* Relativistic GM of each body */
  602. ii = 6*FMASS;
  603. for( i=FMASS; i<NTOTAL; i++ )
  604.     {
  605.     vx = y[ii]; /* velocity */
  606.     s = vx * vx;
  607.     vx = y[ii+2];
  608.     s += vx * vx;
  609.     vx = y[ii+4];
  610.     s += vx * vx;  /* square of velocity */
  611.     for( j=FMASS; j<NTOTAL; j++ )
  612.         {
  613.         if( j == i )
  614.             continue;
  615.         s -= GMs[j]*Rij[i][j];
  616.         }
  617. /*
  618.  *    mustar[i] = GMs[i] * (One + csqi * s);
  619.  */
  620.     mustar[i] = GMs[i] + GMs[i] * csqi * s;
  621.     ii += 6;
  622.     }
  623. /* Compute center of mass with Sun omitted. */
  624. xx = Zero;
  625. yy = Zero;
  626. zz = Zero;
  627. vx = Zero;
  628. vy = Zero;
  629. vz = Zero;
  630. ii = 6*FMASS;
  631. for( i=FMASS; i<NTOTAL; i++ )
  632.     {
  633.     if( i != ISUN )
  634.         {
  635.         s = mustar[i];
  636.         xx += y[ii+1] * s; /* position */
  637.         yy += y[ii+3] * s;
  638.         zz += y[ii+5] * s;
  639.         vx += y[ii] * s; /* velocity */
  640.         vy += y[ii+2] * s;
  641.         vz += y[ii+4] * s;
  642.         }
  643.     ii += 6;
  644.     }
  645. /* Evaluate the Sun so that the center = 0. */
  646. s = mustar[ISUN];
  647. xx = -xx/s;
  648. yy = -yy/s;
  649. zz = -zz/s;
  650. vx = -vx/s;
  651. vy = -vy/s;
  652. vz = -vz/s;
  653. jj = 6*ISUN;
  654. y[jj+1] = xx;
  655. y[jj+3] = yy;
  656. y[jj+5] = zz;
  657. y[jj] = vx;
  658. y[jj+2] = vy;
  659. y[jj+4] = vz;
  660. v[jj+1] = vx;
  661. v[jj+3] = vy;
  662. v[jj+5] = vz;
  663. /* Adjust the table of distances between bodies i and j.
  664.  * Note, most of this work was already done by func().
  665.  * Only the entries involving the Sun need to be changed.
  666.  */
  667. ii = 6*FMASS;
  668. for( j=FMASS; j<NTOTAL; j++ )
  669.     {
  670.     if( j != ISUN )
  671.         {
  672.         vx = xx - y[ii+1];
  673.         s = vx * vx;
  674.         vx = yy - y[ii+3];
  675.         s += vx * vx;
  676.         vx = zz - y[ii+5];
  677.         s += vx * vx;
  678.         vx = SQRT(s);
  679.         vy = One/vx;
  680.         Rij[ISUN][j] = vy;
  681.         Rij[j][ISUN] = vy;
  682.         s = One/(vx*s);
  683.         Rij3[ISUN][j] = s;
  684.         Rij3[j][ISUN] = s;
  685.         }
  686.     ii += 6;
  687.     }
  688.  } /* iteration */
  689.  
  690. }
  691.  
  692.  
  693. #if MOON
  694. /* Convert EMB, M to barycentric Earth and Moon vectors.
  695.  * ymoon[] declared externally.
  696.  */
  697. fromemb( yinp, y )
  698. DOUBLE yinp[], y[];
  699. {
  700. DOUBLE zd, yd;
  701. int i, ii, jj;
  702.  
  703. for( i=0; i<(6*NTOTAL); i++ )
  704.     y[i] = yinp[i];
  705.  
  706. for( i=0; i<6; i++ )
  707.     ymoon[i] = yinp[6*IMOON+i]; /* M */
  708.  
  709. ii = 6*IEARTH;
  710. jj = 6*IMOON;
  711. for( i=0; i<6; i++ )
  712.     {
  713.     zd = yinp[ii+i];        /* EMB */
  714.     yd = yinp[jj+i]/(One+EMRAT);    /* M */
  715.     y[ii+i] = zd -  yd;        /* r_e */
  716.     y[jj+i] = zd + EMRAT * yd;    /* r_m */
  717.     }
  718. }
  719. #endif
  720.  
  721.  
  722.  
  723.  
  724. /* Make table of distances between ith and jth objects
  725.  */
  726. distances(y)
  727. DOUBLE y[];
  728. {
  729. register DOUBLE t, u;
  730. DOUBLE r, xv, yv, zv;
  731. int j, i, jj, ii;
  732.  
  733. jj = 6*(FMASS+1);
  734. for(j=(FMASS+1); j<NTOTAL; j++)
  735.     {
  736.     ii = 6*FMASS;
  737.     xv = y[jj+1];
  738.     yv = y[jj+3];
  739.     zv = y[jj+5];
  740.     for(i=FMASS; i<j; i++)
  741.         {
  742.         t = xv - y[ii+1];
  743.         u = t * t;
  744.         t = yv - y[ii+3];
  745.         u += t * t;
  746.         t = zv - y[ii+5];
  747.         u += t * t;
  748.         r = SQRT(u);
  749.         t = One/r;
  750.         Rij[j][i] = t;
  751.         Rij[i][j]  = t;
  752.         t = One/(r*u);
  753.         Rij3[j][i] = t;
  754.         Rij3[i][j] = t;
  755.         ii += 6;
  756.         }
  757.     jj += 6;
  758.     }
  759. #if MOON
  760. /* Take the input M vector for distance from Earth to Moon.
  761.  * ymoon[] set by previous call to fromemb().
  762.  */
  763. t = ymoon[1];
  764. u = t * t;
  765. t = ymoon[3];
  766. u += t * t;
  767. t = ymoon[5];
  768. u += t * t;
  769. r = SQRT(u);
  770. t = One/r;
  771. Rij[IEARTH][IMOON] = t;
  772. Rij[IMOON][IEARTH] = t;
  773. t = One/(r*u);
  774. Rij3[IEARTH][IMOON] = t;
  775. Rij3[IMOON][IEARTH] = t;
  776. #endif
  777. }
  778.  
  779.  
  780. /* Write date and solution vector to output file.
  781.  */
  782. #if OUTDOUBLE
  783. double yout[6*(ISUN+1)]; /* Output array for disc file. */
  784. #else
  785. DOUBLE yout[6*(ISUN+1)];
  786. #endif
  787.  
  788. wfile( t, y )
  789. DOUBLE t;
  790. DOUBLE y[];
  791. {
  792. DOUBLE p[3], v[3];
  793. #if OUTDOUBLE
  794. double dt;
  795. #else
  796. DOUBLE dt;
  797. #endif
  798. int i, j, k;
  799.  
  800. dt = t;
  801. /*for( i=0; i<6*(ISUN+1); i++ ) */
  802. for( i=0; i<6*FMASS; i++ )
  803.     yout[i] = y[i];
  804.  
  805. /* Rotate to DE200 coordinate system
  806.  */
  807. for( i=FMASS; i<=ISUN; i++ )
  808.     {
  809.     k = 6 * i;
  810.     for( j=0; j<3; j++ )
  811.         {
  812.         v[j] = y[k+j+j];
  813.         p[j] = y[k+j+j+1];
  814.         }
  815. #if DE118
  816.     r118_200(v);
  817.     r118_200(p);
  818. #endif
  819.     for( j=0; j<3; j++ )
  820.         {
  821.         yout[k+j+j] = v[j];
  822.         yout[k+j+j+1] = p[j];
  823.         }
  824.     }
  825.  
  826. #if OUTDOUBLE
  827. write( fd, &dt, sizeof(double) );
  828. write( fd, yout, 6*(ISUN+1)*sizeof(double) );
  829. #else
  830. write( fd, &dt, sizeof(DOUBLE) );
  831. write( fd, yout, 6*(ISUN+1)*sizeof(DOUBLE) );
  832. #endif
  833. }
  834.